Simulation Result

Comparison of Causal Discovery Agorithms

1 Simualtion Design

1.1 Sparse 5 nodes model

See code here.

par(oma=c(0, 0, 4, 0))
## Specify B matrix
B5sparse = matrix(c(0, 0, 0, 0, 0,
                 1, 0, 0.8, 0, 0,
                 0, 0, 0, 0.9, 0,
                 0, 0.7, 0, 0, 1.5,
                 0, 0, 0, 0, 0), 5, 5, byrow = T)

colnames(B5sparse) <- c("X1", "X2", "X3", "X4", "X5")

# specify layout
layout5 = matrix(c(0,1,
                   0,0,
                   1,-1,
                   2,0,
                   2,1),5,2,byrow = T)

true5psparse <- qgraph::qgraph(t(B5sparse), layout=layout5, labels = colnames(B5sparse), theme="colorblind", vsize = 8, asize = 5)
title("True DCG",  font.main = 1, cex.main = 2, line = 2, outer=TRUE)
## True PAG
truepag_5psparse <- matrix(c(
       0, 2, 2, 0, 0,
       3, 0, 3, 3, 3,
       3, 3, 0, 3, 0,
       0, 3, 3, 0, 3,
       0, 2, 0, 2, 0), 5, 5, byrow = TRUE)
dimnames(truepag_5psparse) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(truepag_5psparse)
title(main = "True PAG", font.main = 1, cex.main = 2, line = 2, outer=TRUE)

1.2 Dense 5 nodes model

See code here.

par(oma=c(0, 0, 4, 0))
## Specify B matrix
B5dense = matrix(c(0, 0, 0, 0, 0,
                    1, 0, 0.8, 0, 0,
                    0, 0, 0, 0.9, 0,
                    0, 0.7, 0, 0, 1.5,
                    1, 0, 0, 0, 0), 5, 5, byrow = T)

colnames(B5dense) <- c("X1", "X2", "X3", "X4", "X5")

true5pdense <- qgraph(t(B5dense), layout=layout5, labels = colnames(B5dense), theme="colorblind", vsize = 8, asize = 5)
title("True DCG",  font.main = 1, cex.main = 2, line = 2, outer=TRUE)
## True PAG
truepag_5pdense <- matrix(c(
  0, 2, 2, 0, 2,
  3, 0, 3, 3, 3,
  3, 3, 0, 3, 0,
  0, 3, 3, 0, 3,
  3, 2, 0, 2, 0), 5, 5, byrow = TRUE)
dimnames(truepag_5pdense) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(truepag_5pdense)
title(main = "True PAG", font.main = 1, cex.main = 2, line = 2, outer=TRUE)

1.3 Sparse 10 nodes model

See code here.

par(oma=c(0, 0, 4, 0))
## Specify B matrix
B10sparse = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                  0.4, 0.8, 0, 0, 0, 0, 0, 0, 0, 0, 
                  0, 0, 0.7, 0, 0, 0.9, 0, 0, 0, 0, 
                  0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
                  0, 0, 0, 0, 0.2, 0, 0.5, 0, 0, 0, 
                  0, 0, 0, 0, 0, 0, 0, 1, 0.8, 0, 
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0.4, 
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 10, 10, byrow = T)

dimnames(B10sparse) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))

# specify layout
layout10 = matrix(c(0,1,
                      2,1,
                      1,0,
                      2,-1,
                      3,0,
                      4, -1,
                      5, 0,
                      6, -1,
                      4, 1,
                      7, 1),10,2,byrow = T)

true10psparse <- qgraph(t(B10sparse), layout = layout10, labels = colnames(B10sparse), theme="colorblind", vsize = 8, asize = 5)
title("True DCG",  font.main = 1, cex.main = 2, line = 2, outer=TRUE)
## True PAG
truepag_10psparse <- matrix(c(
         0, 0, 2, 0, 0, 0, 0, 0, 0, 0,
         0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 
         3, 3, 0, 2, 0, 0, 0, 0, 0, 0, 
         0, 0, 3, 0, 3, 3, 0, 0, 0, 0, 
         0, 0, 0, 3, 0, 3, 0, 0, 0, 0, 
         0, 0, 0, 3, 3, 0, 3, 0, 0, 0, 
         0, 0, 0, 0, 0, 2, 0, 3, 3, 0, 
         0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 
         0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 
         0, 0, 0, 0, 0, 0, 0, 0, 2, 0), 10, 10, byrow = T)

dimnames(truepag_10psparse) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(truepag_10psparse)
title(main = "True PAG", font.main = 1, cex.main = 2, line = 2, outer=TRUE)

1.4 Dense 10 nodes model

See code here.

par(oma=c(0, 0, 4, 0))
## Specify B matrix
B10dense = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                     0.4, 0.8, 0, 0, 0, 0, 0, 0, 0, 0, 
                     0, 0, 0.7, 0, 0, 0.9, 0, 0, 0, 0, 
                     0, 0.4, 0, 1, 0, 0, 0, 0, 0, 0, 
                     0, 0, 0, 0, 0.9, 0, 0.5, 0, 0, 0, 
                     0, 0, 0, 0, 0, 0, 0, 1, 0.8, 1, 
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6, 
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0.4, 
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 10, 10, byrow = T)

colnames(B10dense) <- c(paste("X", 1:10, sep=""))

# specify layout
layout10 = matrix(c(0,1,
                    2,1,
                    1,0,
                    2,-1,
                    3,0,
                    4, -1,
                    5, 0,
                    6, -1,
                    4, 1,
                    7, 1),10,2,byrow = T)

true10pdense <- qgraph(t(B10dense), layout = layout10, labels = colnames(B10dense), theme="colorblind", vsize = 8, asize = 5)
title("True DCG",  font.main = 1, cex.main = 2, line = 2, outer=TRUE)
## True PAG
truepag_10pdense <- matrix(c(
  0, 0, 2, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 2, 2, 2, 0, 0, 0, 0, 0, 
  3, 3, 0, 2, 0, 2, 0, 0, 0, 0, 
  0, 3, 3, 0, 3, 3, 0, 0, 0, 0, 
  0, 3, 0, 3, 0, 3, 3, 0, 0, 0, 
  0, 0, 0, 3, 3, 0, 3, 0, 0, 0, 
  0, 0, 0, 0, 0, 2, 0, 3, 3, 3, 
  0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 
  0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 
  0, 0, 0, 0, 0, 0, 2, 2, 2, 0), 10, 10, byrow = T)

dimnames(truepag_10pdense) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(truepag_10pdense)
title(main = "True PAG", font.main = 1, cex.main = 2, line = 2, outer=TRUE)

1.5 5 nodes with a LV sparse model

See code here.

par(oma=c(0, 0, 4, 0))
## Specify B matrix
B5_lv = matrix(c(0, 0, 0, 0, 0, 1,
                 0, 0, 0.4, 0, 0, 1,
                 0, 0, 0, 0.5, 0,0,
                 0, 0.7, 0, 0, 1.5,0,
                 0, 0, 0, 0, 0,0,
                 0,0,0,0,0,0), 6, 6, byrow = T)

colnames(B5_lv) <- c("X1", "X2", "X3", "X4", "X5", "L1")
# specify layout
layout5_lv = matrix(c(0,1,
                      0,0,
                      1,-1,
                      2,0,
                      2,1,
                      -1, 0.5),6,2,byrow = T)

true5p_lv <- qgraph(t(B5_lv), layout=layout5_lv, labels = colnames(B5_lv), theme="colorblind")
title("True DCG",  font.main = 1, cex.main = 2, line = 2, outer=TRUE)
## True PAG
truepag_5pLV <- matrix(c(
  0, 2, 2, 0, 0,
  2, 0, 3, 3, 3,
  2, 3, 0, 3, 0,
  0, 3, 3, 0, 3,
  0, 2, 0, 2, 0), 5, 5, byrow = TRUE)

dimnames(truepag_5pLV) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(truepag_5pLV)
title(main = "True PAG", font.main = 1, cex.main = 2, line = 2, outer=TRUE)

1.6 5 nodes with a LV dense model

See code here.

par(oma=c(0, 0, 4, 0))
## Specify B matrix
B5_lvdense = matrix(c(0, 0, 0, 0, 0, 1,
                       0, 0, 0.4, 0, 0, 1,
                       0, 0, 0, 0.5, 0,0,
                       0, 0.7, 0, 0, 1.5,0,
                       0.6, 0, 0, 0, 0,0,
                       0,0,0,0,0,0), 6, 6, byrow = T)

colnames(B5_lvdense) <- c("X1", "X2", "X3", "X4", "X5", "L1")
# specify layout
layout5_lv = matrix(c(0,1,
                      0,0,
                      1,-1,
                      2,0,
                      2,1,
                      -1, 0.5),6,2,byrow = T)

true5p_lvdense <- qgraph(t(B5_lvdense), layout=layout5_lv, labels = colnames(B5_lvdense), theme="colorblind")

title("True DCG",  font.main = 1, cex.main = 2, line = 2, outer=TRUE)
## True PAG
truepag_5pLVdense <- matrix(c(
  0, 2, 2, 0, 2,
  2, 0, 3, 3, 3,
  2, 3, 0, 3, 0,
  0, 3, 3, 0, 3,
  3, 2, 0, 2, 0), 5, 5, byrow = TRUE)

dimnames(truepag_5pLVdense) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(truepag_5pLVdense)
title(main = "True PAG", font.main = 1, cex.main = 2, line = 2, outer=TRUE)

1.7 10 nodes with a LV sparse model

See code here.

par(oma=c(0, 0, 4, 0))
## Specify B matrix
B10_lv = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                  0.4, 0.8, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                  0, 0, 0.7, 0, 0, 0.9, 0, 0, 0, 0, 0,
                  0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
                  0, 0, 0, 0, 0.2, 0, 0.5, 0, 0, 0, 0,
                  0, 0, 0, 0, 0, 0, 0, 1, 0.8, 0, 0,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.8,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0.4, 0,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 11, 11, byrow = T)

colnames(B10_lv) <- c(paste("X", 1:10, sep=""), "L1")

# specify layout
layout10LV = matrix(c(0,1,
                      2,1,
                      1,0,
                      2,-1,
                      3,0,
                      4, -1,
                      5, 0,
                      6, -1,
                      4, 1,
                      7, 1,
                      8, 0),11,2,byrow = T)

true10pLV <- qgraph(t(B10_lv), layout = layout10LV, labels = colnames(B10_lv), theme="colorblind", vsize = 8, asize = 5)
title("True DCG",  font.main = 1, cex.main = 2, line = 2, outer=TRUE)
## True PAG
truepag_10pLV <- matrix(c(
  0, 0, 2, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 
  3, 3, 0, 2, 0, 2, 0, 0, 0, 0, 
  0, 0, 3, 0, 3, 3, 0, 0, 0, 0, 
  0, 0, 0, 3, 0, 3, 3, 0, 0, 0, 
  0, 0, 3, 3, 3, 0, 3, 0, 0, 0, 
  0, 0, 0, 0, 2, 2, 0, 3, 3, 0, 
  0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 
  0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 
  0, 0, 0, 0, 0, 0, 0, 2, 2, 0), 10, 10, byrow = T)

dimnames(truepag_10pLV) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(truepag_10pLV)
title(main = "True PAG", font.main = 1, cex.main = 2, line = 2, outer=TRUE)

1.8 10 nodes with a LV dense model

See code here.

par(oma=c(0, 0, 4, 0))
## Specify B matrix
B10_lvdense = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                  0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                  0.4, 0.8, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                  0, 0, 0.7, 0, 0, 0.9, 0, 0, 0, 0, 0,
                  0, 0, 0.6, 1, 0, 0, 0, 0, 0, 0, 0,
                  0, 0, 0, 0, 0.2, 0, 0.5, 0, 0, 0, 0,
                  0, 0, 0, 0, 0, 0, 0, 1, 0.8, 0.6, 0,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.8,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0.4, 0,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2,
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 11, 11, byrow = T)

colnames(B10_lvdense) <- c(paste("X", 1:10, sep=""), "L1")

# specify layout
layout10LV = matrix(c(0,1,
                      2,1,
                      1,0,
                      2,-1,
                      3,0,
                      4, -1,
                      5, 0,
                      6, -1,
                      4, 1,
                      7, 1,
                      8, 0),11,2,byrow = T)

true10pLVdense <- qgraph(t(B10_lvdense), layout = layout10LV, labels = colnames(B10_lvdense), theme="colorblind")

title("True DCG",  font.main = 1, cex.main = 2, line = 2, outer=TRUE)
## True PAG
truepag_10pLVdense <- matrix(c(
  0, 2, 2, 0, 0, 0, 0, 0, 0, 0,
  3, 0, 2, 0, 0, 0, 0, 0, 0, 0, 
  3, 3, 0, 2, 2, 2, 0, 0, 0, 0, 
  0, 0, 3, 0, 3, 3, 0, 0, 0, 0, 
  0, 0, 3, 3, 0, 3, 3, 0, 0, 0, 
  0, 0, 3, 3, 3, 0, 3, 0, 0, 0, 
  0, 0, 0, 0, 2, 2, 0, 3, 3, 3, 
  0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 
  0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 
  0, 0, 0, 0, 0, 0, 2, 2, 2, 0), 10, 10, byrow = T)

dimnames(truepag_10pLVdense) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(truepag_10pLVdense)
title(main = "True PAG", font.main = 1, cex.main = 2, line = 2, outer=TRUE)

---

2 Simulation Result

See code here.

## Precision & Recall
# compute the average precision and sd
results <- list(res_ccd5psparse, res_fci5psparse, res_cci5psparse, res_ccd10psparse, res_fci10psparse, res_cci10psparse, res_ccd5pdense, res_fci5pdense, res_cci5pdense, res_ccd10pdense, res_fci10pdense, res_cci10pdense, res_ccd5pLVsparse, res_fci5pLVsparse, res_cci5pLVsparse, res_ccd5pLVdense, res_fci5pLVdense, res_cci5pLVdense, res_ccd10pLVsparse,  res_fci10pLVsparse, res_cci10pLVsparse, res_ccd10pdense,  res_fci10pLVdense, res_cci10pLVdense) %>% 
  # transpose df
  map(~ sjmisc::rotate_df(.x) %>%
  # add sample size (N) info
  rename_with(~paste0(.x, "N = ", rep(N, each=8)))  %>%
  # think about how to deal with NAs or do I want to define sth. else instead of NAs.
  #na.omit(.x) %>% 
  summarise(across(everything(.), list(mean = ~mean(., na.rm=T), sd = ~sd(., na.rm=T))))) %>% 
  bind_rows() %>% 
  mutate(algorithm = rep(c("ccd", "fci", "cci"), 8),
         condition = rep(c("5p_sparse", "10p_sparse", "5p_dense", "10p_dense", "5p_LVsparse", "5p_LVdense", "10p_LVsparse", "10p_LVdense"), each=3)) %>%
  # brings the algorithm and condition names first
  relocate(where(is.character), .before = where(is.numeric)) %>% 
  # convert itto a long format
  tidyr::pivot_longer(!c(algorithm, condition), names_to = "metric", values_to = "value") %>% 
  # Add sample size column (N) & clean up the column name
  mutate(N = stringr::str_extract(metric, "(?<=[N =])\\d+"),
         metric = stringr::str_replace_all(metric, "[0-9.]+|[N =]", "")) 

## Uncertainty
uncertainties <- bind_rows("ccd_5p-sparse" = uncer_ccd5psparse, "fci_5p-sparse" = uncer_fci5psparse, "cci_5p-sparse"=uncer_cci5psparse, "ccd_10p-sparse"=uncer_ccd10psparse, "fci_10p-sparse" = uncer_fci10psparse, "cci_10p-sparse" = uncer_cci10psparse, "ccd_5p-dense"=uncer_ccd5pdense, "fci_5p-dense"=uncer_fci5pdense, "cci_5p-dense"=uncer_cci5pdense, "ccd_10p-dense"=uncer_ccd10pdense, "fci_10p-dense"=uncer_fci10pdense, "cci_10p-dense"=uncer_cci10pdense, "ccd_5p-LVsparse"=uncer_ccd5pLVsparse, "fci_5p-LVsparse"=uncer_fci5pLVsparse, "cci_5p-LVsparse"=uncer_cci5pLVsparse, "ccd_10p-LVsparse"=uncer_ccd10pLVsparse, "fci_10p-LVsparse"=uncer_fci10pLVsparse, "cci_10p-LVsparse"=uncer_cci10pLVsparse,
"ccd_5p-LVdense"=uncer_ccd5pLVdense, "fci_5p-LVdense"=uncer_fci5pLVdense, "cci_5p-LVdense"=uncer_cci5pLVdense, "ccd_10p-LVdense"=uncer_ccd10pLVdense, "fci_10p-LVdense"=uncer_fci10pLVdense, "cci_10p-LVdense"=uncer_cci10pLVdense,.id="id") %>% 
  group_by(id) %>% 
  summarise_all(list(mean = mean, sd = sd)) %>%  
  mutate(algorithm = stringr::str_split(id, "_", simplify = T)[,1],
         condition = stringr::str_split(id, "_", simplify = T)[,2]) %>% 
  tidyr::pivot_longer(!c(algorithm, condition, id), names_to = "name", values_to = "value") %>% 
  mutate(N = stringr::str_extract(stringr::str_split(name, "_", simplify = T)[,1], "(\\d)+"),
         statistics = stringr::str_split(name, "_", simplify = T)[,2]) %>% 
  dplyr::select(-id, -name) %>%  relocate(where(is.character), .before = where(is.numeric))



## SHD
SHDs <- bind_rows("ccd_5p-sparse" = SHD_ccd5psparse, "fci_5p-sparse" = SHD_fci5psparse, "cci_5p-sparse"=SHD_cci5psparse, "ccd_10p-sparse"= SHD_ccd10psparse, "fci_10p-sparse" = SHD_fci10psparse, "cci_10p-sparse" = SHD_cci10psparse, "ccd_5p-dense"= SHD_ccd5pdense, "fci_5p-dense"=SHD_fci5pdense, "cci_5p-dense"=SHD_cci5pdense, "ccd_10p-dense"= SHD_ccd10pdense, "fci_10p-dense"=SHD_fci10pdense, "cci_10p-dense"=SHD_cci10pdense, "ccd_5p-LVsparse"=SHD_ccd5pLVsparse, "fci_5p-LVsparse"=SHD_fci5pLVsparse, "cci_5p-LVsparse"=SHD_cci5pLVsparse, "ccd_10p-LVsparse"=SHD_ccd10pLVsparse, "fci_10p-LVsparse"=SHD_fci10pLVsparse, "cci_10p-LVsparse"=SHD_cci10pLVsparse, 
"ccd_5p-LVdense"=SHD_ccd5pLVdense, "fci_5p-LVdense"=SHD_fci5pLVdense, "cci_5p-LVdense"=SHD_cci5pLVdense, "ccd_10p-LVdense"=SHD_ccd10pLVdense, "fci_10p-LVdense"=SHD_fci10pLVdense, "cci_10p-LVdense"=SHD_cci10pLVdense, .id="id") %>% 
  group_by(id) %>% 
  summarise_all(list(mean = mean, sd = sd)) %>%  
  mutate(algorithm = stringr::str_split(id, "_", simplify = T)[,1],
         condition = stringr::str_split(id, "_", simplify = T)[,2]) %>% 
  tidyr::pivot_longer(!c(algorithm, condition, id), names_to = "name", values_to = "value") %>% 
  mutate(N = stringr::str_extract(stringr::str_split(name, "_", simplify = T)[,1], "(\\d)+"),
         statistics = stringr::str_split(name, "_", simplify = T)[,2]) %>% 
  dplyr::select(-id, -name) %>%  relocate(where(is.character), .before = where(is.numeric)) 

See code here.

## Specify my custom theme
MyTheme <-  theme(plot.title = element_blank(),
                  plot.subtitle = element_text( face = "italic"),
                  axis.text=element_text(face = "bold"),
                  legend.text = element_text(face = "bold"))

## ========================
## precision plots
## ========================
precision_plots <- c(unique(results$condition)) %>% 
  map(~
results %>% 
  filter(condition == .x & grepl("average_precision", metric)) %>% 
  tidyr::pivot_wider(names_from = metric, values_from=value) %>% 
  ggplot(aes(x= factor(N, levels = c("50", "150", "500", "1000", "5000")), y=average_precision_mean, group = algorithm, colour = algorithm, fill=algorithm)) +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  #geom_errorbar(aes(ymin=average_precision_mean-qnorm(0.975)*average_precision_sd/sqrt(as.numeric(N)), ymax=average_precision_mean+qnorm(0.975)*average_precision_sd/sqrt(as.numeric(N))), width=0.1) +
  geom_ribbon(aes(ymin=average_precision_mean-qnorm(0.975)*average_precision_sd/sqrt(as.numeric(N)), ymax=average_precision_mean+qnorm(0.975)*average_precision_sd/sqrt(as.numeric(N))), alpha=0.2,lwd=0) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  labs(x="N", y="", title = "", subtitle = .x) +
  theme_classic() + MyTheme
)

## ========================
## recall plots
## ========================
recall_plots <- c(unique(results$condition)) %>% 
  map(~
results %>% 
  filter(condition == .x & grepl("average_recall", metric)) %>% 
  tidyr::pivot_wider(names_from = metric, values_from=value) %>% 
  ggplot(aes(x= factor(N, levels = c("50", "150", "500", "1000", "5000")), y=average_recall_mean, group = algorithm, colour = algorithm, fill= algorithm)) +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  #geom_errorbar(aes(ymin=average_recall_mean-qnorm(0.975)*average_recall_sd/sqrt(as.numeric(N)), ymax=average_recall_mean+qnorm(0.975)*average_recall_sd/sqrt(as.numeric(N))), width=0.1) +
  geom_ribbon(aes(ymin=average_recall_mean-qnorm(0.975)*average_recall_sd/sqrt(as.numeric(N)), ymax=average_recall_mean+qnorm(0.975)*average_recall_sd/sqrt(as.numeric(N))), alpha=0.2,lwd=0) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  labs(x="N", y="", title = "", subtitle = .x) +
  theme_classic() + MyTheme
)

## ========================
## uncertainty plots
## ========================
## uncertainty (uncertainty rate of fci and cci are exactly the same!)
uncertainty_plots <- c(unique(uncertainties$condition)) %>% 
  map(~
uncertainties %>%
  filter(condition == .x) %>% 
  tidyr::pivot_wider(names_from = statistics, values_from=value) %>% 
  ggplot(aes(x= factor(N, levels = c("50", "150", "500", "1000", "5000")), y=mean, group = algorithm, colour = algorithm, fill=algorithm)) +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  #geom_errorbar(aes(ymin=mean-qnorm(0.975)*sd/sqrt(as.numeric(N)), ymax=mean+qnorm(0.975)*sd/sqrt(as.numeric(N))), width=0.1) +
  geom_ribbon(aes(ymin=mean-qnorm(0.975)*sd/sqrt(as.numeric(N)), ymax=mean+qnorm(0.975)*sd/sqrt(as.numeric(N))), alpha=0.2,lwd=0) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  labs(x="N", y="", title = "", subtitle = .x) +
  theme_classic() + MyTheme
)


## ========================
## SHD plots 
## ========================
SHD_plots <- c(unique(SHDs$condition)) %>% 
  map(~ 
SHDs %>%
  filter(condition == .x) %>% 
  tidyr::pivot_wider(names_from = statistics, values_from=value) %>% 
  ggplot(aes(x= factor(N, levels = c("50", "150", "500", "1000", "5000")), y=mean, group = algorithm, colour = algorithm, fill = algorithm)) +
  geom_line(aes(group = algorithm)) +
  geom_point() + 
  #geom_errorbar(aes(ymin=mean-qnorm(0.975)*sd/sqrt(as.numeric(N)), ymax=mean+qnorm(0.975)*sd/sqrt(as.numeric(N))), width=0.1) +
  geom_ribbon(aes(ymin=mean-qnorm(0.975)*sd/sqrt(as.numeric(N)), ymax=mean+qnorm(0.975)*sd/sqrt(as.numeric(N))), alpha=0.2,lwd=0) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  labs(x="N", y="", title = "", subtitle = .x) +
  theme_classic() + MyTheme
)


SHDs |> 
  filter(condition == "10p-LVdense") |>
  tidyr::pivot_wider(names_from = statistics, values_from = value)

# A tibble: 15 × 5
   algorithm condition   N      mean    sd
   <chr>     <chr>       <chr> <dbl> <dbl>
 1 ccd       10p-LVdense 50     31.6  2.88
 2 ccd       10p-LVdense 150    27.6  2.56
 3 ccd       10p-LVdense 500    25.5  3.33
 4 ccd       10p-LVdense 1000   22.3  4.11
 5 ccd       10p-LVdense 5000   17.2  4.19
 6 cci       10p-LVdense 50     34.9  2.02
 7 cci       10p-LVdense 150    29.9  1.94
 8 cci       10p-LVdense 500    26.6  1.74
 9 cci       10p-LVdense 1000   25.5  2.24
10 cci       10p-LVdense 5000   19.4  3.24
11 fci       10p-LVdense 50     36.9  2.48
12 fci       10p-LVdense 150    30.3  2.62
13 fci       10p-LVdense 500    26.7  2.02
14 fci       10p-LVdense 1000   25.6  2.13
15 fci       10p-LVdense 5000   21.2  1.77

See code here.

# combine plots
# precision plot
ggarrange(plotlist = precision_plots,
                    ncol = 2, nrow = 4, common.legend = TRUE, legend = "bottom") %>%
  annotate_figure(top = text_grob("Precision", face = "bold", size = 14, family = "Palatino"))
# recall plot
ggarrange(plotlist = recall_plots,
                    ncol = 2, nrow = 4, common.legend = TRUE, legend = "bottom") %>%
  annotate_figure(top = text_grob("Recall", face = "bold", size = 14, family = "Palatino"))
# uncertainty plot
ggarrange(plotlist = uncertainty_plots,
                    ncol = 2, nrow = 4, common.legend = TRUE, legend = "bottom") %>%
  annotate_figure(top = text_grob("Uncertainty", face = "bold", size = 14, family = "Palatino"))
# shd plot
ggarrange(plotlist = SHD_plots,
                    ncol = 2, nrow = 4, common.legend = TRUE, legend = "bottom") %>% 
annotate_figure(top = text_grob("SHD", face = "bold", size = 14, family = "Palatino"))

Figure 1. Simuation results_varying sample sizes

Figure 2. Simuation results_varying sample sizes

Figure 3. Simuation results_varying sample sizes

Figure 4. Simuation results_varying sample sizes

See code here.

## ========================
## WITHOUT LV CONDITION
## ========================

## WITHOUT LV CONDITION: PRECISION
p1 <- results %>% 
  # exclude LV conditions
  filter(!grepl("LV", condition), N == "1000") %>% 
  tidyr::pivot_wider(names_from = metric, values_from=value) %>% 
ggplot(aes(x= factor(condition, levels = c("5p_sparse", "5p_dense", "10p_sparse", "10p_dense")), y=average_precision_mean, group = algorithm, colour = algorithm)) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  # in this case, N = 1000 (hence, sqrt(1000))
  geom_errorbar(aes(ymin=average_precision_mean-qnorm(0.975)*average_precision_sd/sqrt(as.numeric(N)), ymax=average_precision_mean+qnorm(0.975)*average_precision_sd/sqrt(as.numeric(N))), width=0.1) +
  labs(x="", y="",title = "Without a Latent Variable Condition", subtitle = "Without LV_PRECISION") +
  theme_classic() + MyTheme

## WITHOUT LV CONDITION: RECALL
p2 <- results %>% 
  # exclude LV conditions
  filter(!grepl("LV", condition), N == "1000") %>% 
  tidyr::pivot_wider(names_from = metric, values_from=value) %>% 
  ggplot(aes(x= factor(condition, levels = c("5p_sparse", "5p_dense", "10p_sparse", "10p_dense")), y=average_recall_mean, group = algorithm, colour = algorithm)) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  # in this case, N = 1000 (hence, sqrt(1000))
  geom_errorbar(aes(ymin=average_recall_mean-qnorm(0.975)*average_recall_sd/sqrt(as.numeric(N)), ymax=average_recall_mean+qnorm(0.975)*average_recall_sd/sqrt(as.numeric(N))), width=0.1) +
  labs(x="", y="", title = "", subtitle = "Without LV_RECALL") +
  theme_classic() + MyTheme
  

## WITHOUT LV CONDITION: UNCERTAINTY
p3 <- uncertainties %>% 
  # exclude LV conditions
  filter(!grepl("LV", condition), N == "1000") %>% 
  tidyr::pivot_wider(names_from = statistics, values_from=value) %>% 
  ggplot(aes(x= factor(condition, levels = c("5p-sparse", "5p-dense", "10p-sparse", "10p-dense")), y=mean, group = algorithm, colour = algorithm)) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  # in this case, N = 1000 (hence, sqrt(1000))
  geom_errorbar(aes(ymin=mean-qnorm(0.975)*sd/sqrt(1000), ymax=mean+qnorm(0.975)*sd/sqrt(1000)), width=0.1) +
  labs(x="", y="", title = "", subtitle = "Without LV_UNCERTAINTY") +
  theme_classic() + MyTheme
  

## WITHOUT LV CONDITION: SHD
p4 <- SHDs %>% 
  # exclude LV conditions
  filter(!grepl("LV", condition), N == "1000") %>% 
  tidyr::pivot_wider(names_from = statistics, values_from=value) %>% 
  ggplot(aes(x= factor(condition, levels = c("5p-sparse", "5p-dense", "10p-sparse", "10p-dense")), y=mean, group = algorithm, colour = algorithm)) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  # in this case, N = 1000 (hence, sqrt(1000)) 
  geom_errorbar(aes(ymin=mean-qnorm(0.975)*sd/sqrt(1000), ymax=mean+qnorm(0.975)*sd/sqrt(1000)), width=0.1) +
  labs(x="", y="", title = "", subtitle = "Without LV_SHD") +
  theme_classic() + MyTheme

  


## ========================
## WITH LV CONDITION
## ========================

## WITH LV CONDITION: PRECISION
p5 <- results %>% 
  # exclude LV conditions
  filter(grepl("LV", condition), N == "1000") %>% 
  tidyr::pivot_wider(names_from = metric, values_from=value) %>% 
  ggplot(aes(x= factor(condition, levels = c("5p_LVsparse", "10p_LVsparse", "5p_LVdense", "10p_LVdense")), y=average_precision_mean, group = algorithm, colour = algorithm)) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  # in this case, N = 1000 (hence, sqrt(1000))
  geom_errorbar(aes(ymin=average_precision_mean-qnorm(0.975)*average_precision_sd/sqrt(as.numeric(N)), ymax=average_precision_mean+qnorm(0.975)*average_precision_sd/sqrt(as.numeric(N))), width=0.1) +
  labs(x="", y="", title = "With a Latent Variable Condition", subtitle = "With LV_PRECISION") +
  theme_classic() + MyTheme


## WITH LV CONDITION: RECALL
p6 <- results %>% 
  # exclude LV conditions
  filter(grepl("LV", condition), N == "1000") %>% 
  tidyr::pivot_wider(names_from = metric, values_from=value) %>% 
  ggplot(aes(x= factor(condition, levels = c("5p_LVsparse", "10p_LVsparse", "5p_LVdense", "10p_LVdense")), y=average_recall_mean, group = algorithm, colour = algorithm)) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  # in this case, N = 1000 (hence, sqrt(1000))
  geom_errorbar(aes(ymin=average_recall_mean-qnorm(0.975)*average_recall_sd/sqrt(as.numeric(N)), ymax=average_recall_mean+qnorm(0.975)*average_recall_sd/sqrt(as.numeric(N))), width=0.1) +
  labs(x="", y="", title = "", subtitle = "With LV_RECALL") +
  theme_classic() + MyTheme
  

## WITH LV CONDITION: UNCERTAINTY
p7 <- uncertainties %>% 
  # exclude LV conditions
  filter(grepl("LV", condition), N == "1000") %>% 
  tidyr::pivot_wider(names_from = statistics, values_from=value) %>% 
  ggplot(aes(x= factor(condition, levels = c("5p-LVsparse", "10p-LVsparse", "5p-LVdense", "10p-LVdense")), y=mean, group = algorithm, colour = algorithm)) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  # in this case, N = 1000 (hence, sqrt(1000))
  geom_errorbar(aes(ymin=mean-qnorm(0.975)*sd/sqrt(as.numeric(N)), ymax=mean+qnorm(0.975)*sd/sqrt(as.numeric(N))), width=0.1) +
  labs(x="", y="", title = "", subtitle = "With LV_UNCERTAINTY") +
  theme_classic() + MyTheme
  

## WITH LV CONDITION: SHD
p8 <- SHDs %>% 
  # exclude LV conditions
  filter(grepl("LV", condition), N == "1000") %>% 
  tidyr::pivot_wider(names_from = statistics, values_from=value) %>% 
  ggplot(aes(x= factor(condition, levels = c("5p-LVsparse", "10p-LVsparse", "5p-LVdense", "10p-LVdense")), y=mean, group = algorithm, colour = algorithm)) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  geom_line(aes(group = algorithm)) +
  geom_point() +
  # in this case, N = 1000 (hence, sqrt(1000))
  geom_errorbar(aes(ymin=mean-qnorm(0.975)*sd/sqrt(as.numeric(N)), ymax=mean+qnorm(0.975)*sd/sqrt(as.numeric(N))), width=0.1) +
  labs(x="", y="", title = "", subtitle = "With LV_SHD") +
  theme_classic() + MyTheme
  

# combine plots
ggarrange(p1, p5, p2, p6, p3, p7, p4, p8,
                    ncol = 2, nrow = 4, common.legend = TRUE, legend = "bottom")

Figure 5. Simuation results

See code here.

## plot the results

times %>%
  ggplot(aes(x=factor(condition, levels= c("5psparse", "5pdense", "10psparse", "10pdense", "5pLVsparse","5pLVdense", "10pLVsparse","10pLVdense")), y = log(time), col= factor(algorithm))) +
  geom_boxplot(position = "dodge",   outlier.size = 0.8, outlier.alpha = 0.2) + theme_classic() +
  # scale_x_discrete(name ="Condition",
  #                  labels=c("", "5p-sparse", "", "","5p-dense","","", "10p-sparse","","","10p-dense","","","5p-LV","","","10p-LV","")) +
  scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
  labs(y = " log(ms)", x = "conditions", title = "Algorithm Running Time", subtitle = "Time in milliseconds (ms)") +
  theme(axis.text.x = element_text(face = "bold", margin = margin(t = 13), size=10),
        legend.position="bottom",
        legend.text = element_text(face = "bold"))

Figure 6. Algorithm running time

---

3 Empirical Example

See code here.

# empirical data 408 rows by 26 columns (p = 26)
mcnally <- read.csv("../data/McNally.csv") 
# check the data
#glimpse(mcnally)
skimr::skim(mcnally)

Data summary

Name	mcnally
Number of rows	408
Number of columns	26
_______________________
Column type frequency:
numeric	26
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
onset	0	1	1.20	1.07	0	0	1	2.00	3	▇▆▁▆▃
middle	0	1	1.44	1.07	0	0	2	2.00	3	▆▃▁▇▃
late	0	1	0.81	1.07	0	0	0	2.00	3	▇▂▁▂▂
hypersom	0	1	1.01	0.99	0	0	1	2.00	3	▇▇▁▅▂
sad	0	1	1.55	0.94	0	1	2	2.00	3	▃▇▁▇▅
decappetite	0	1	0.49	0.73	0	0	0	1.00	3	▇▃▁▁▁
incappetite	0	1	0.44	0.87	0	0	0	0.00	3	▇▁▁▁▁
weightloss	0	1	0.50	0.94	0	0	0	1.00	3	▇▁▁▁▁
weightgain	0	1	0.67	1.04	0	0	0	1.00	3	▇▂▁▂▁
concen	0	1	1.48	0.87	0	1	1	2.00	3	▃▇▁▇▂
guilt	0	1	1.56	1.17	0	1	1	3.00	3	▆▇▁▃▇
suicide	0	1	0.63	0.82	0	0	0	1.00	3	▇▅▁▂▁
anhedonia	0	1	1.27	1.05	0	0	1	2.00	3	▇▇▁▅▅
fatigue	0	1	1.33	0.95	0	1	1	2.00	3	▅▇▁▆▃
retard	0	1	0.66	0.81	0	0	0	1.00	3	▇▅▁▂▁
agitation	0	1	1.10	0.93	0	0	1	2.00	3	▅▇▁▃▂
obtime	0	1	2.95	0.89	0	2	3	4.00	4	▁▁▆▇▇
obinterfer	0	1	2.69	0.83	0	2	3	3.00	4	▁▁▆▇▃
obdistress	0	1	2.81	0.76	0	2	3	3.00	4	▁▁▅▇▃
obresist	0	1	1.98	0.93	0	1	2	2.25	4	▁▃▇▃▁
obcontrol	0	1	2.67	0.76	0	2	3	3.00	4	▁▁▆▇▂
comptime	0	1	2.60	0.93	0	2	3	3.00	4	▁▂▇▇▅
compinterf	0	1	2.58	0.88	0	2	3	3.00	4	▁▂▆▇▂
compdis	0	1	2.71	0.84	0	2	3	3.00	4	▁▁▅▇▂
compresis	0	1	2.16	0.89	0	2	2	3.00	4	▁▂▇▅▁
compcont	0	1	2.60	0.76	0	2	3	3.00	4	▁▁▆▇▂

See code here.

# separate dep / ocd symptoms
depression <- mcnally[,1:16]
ocd <- mcnally[,17:26]

See code here.

## run ccd algorithm on entire mcnally (discrete)
ccd_mcnally <- ccdKP(df=mcnally, dataType = "discrete") 
matccd_mcnally <- CreateAdjMat(ccd_mcnally, p = 26)

## run fci algorithm on entire mcnally (discrete)
fci_mcnally <- tetradrunner(algoId = 'fci', df = mcnally,
                          dataType = 'discrete')
matfci_mcnally <- CreateAdjMat(fci_mcnally, p = ncol(mcnally))

## run ccd on depression symptoms
ccd_mcnally_dep <- ccdKP(df=depression, dataType = "discrete") 
mat_mcnally_dep <- CreateAdjMat(ccd_mcnally_dep, p = ncol(depression))

## run fci on depression symptoms
fci_mcdep <- tetradrunner(algoId = 'fci', df = depression,
                          dataType = 'discrete')
matfci_mcdep <- CreateAdjMat(fci_mcdep, p = ncol(depression))

## run ccd on ocd symptoms
ccd_mcnally_ocd <- ccdKP(df=ocd, dataType = "discrete") 
mat_mcnally_ocd <- CreateAdjMat(ccd_mcnally_ocd, p = ncol(ocd))

## run fci on ocd symptoms
fci_mcocd <- tetradrunner(algoId = 'fci', df = ocd,
                          dataType = 'discrete')
matfci_mcocd <- CreateAdjMat(fci_mcocd, p = ncol(ocd))

See code here.

## CCD PAG for the entire Mcnally
plotPAG(ccd_mcnally, matccd_mcnally) 

Figure 7. CCD PAG on entire McNally data

See code here.

## FCI PAG for the entire Mcnally
plotPAG(fci_mcnally, matfci_mcnally) 

Figure 8. FCI PAG on entire McNally data

See code here.

## CCD PAG for the depression symptoms
plotPAG(ccd_mcnally_dep, mat_mcnally_dep) 

## FCI PAG for the depression symptoms
plotPAG(fci_mcdep, matfci_mcdep)

## CCD PAG for the ocd symptoms
plotPAG(ccd_mcnally_ocd, mat_mcnally_ocd) 

## FCI PAG for the ocd symptoms
plotPAG(fci_mcocd, matfci_mcocd)

(a) CCD PAG_depression symptoms

(b) FCI PAG_depression symptoms

(c) CCD PAG_OCD symptoms

(d) FCI PAG_OCD symptoms

Figure 9. Separate PAGs on depression and OCD symptoms

4 Session info

sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.0

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] microbenchmark_1.4.9 CCI.KP_0.1.0         MASS_7.3-58.2       
 [4] furrr_0.3.1          future_1.31.0        magrittr_2.0.3      
 [7] Rgraphviz_2.40.0     graph_1.74.0         BiocGenerics_0.42.0 
[10] DOT_0.1              rcausal_1.2.1        devtools_2.4.5      
[13] usethis_2.1.6        rJava_1.0-6          dplyr_1.1.0         
[16] purrr_1.0.1          ggpubr_0.5.0         ggplot2_3.4.0       
[19] pcalg_2.7-8          qgraph_1.9.3        

loaded via a namespace (and not attached):
  [1] backports_1.4.1     Hmisc_4.7-2         plyr_1.8.8         
  [4] igraph_1.3.5        repr_1.1.6          splines_4.2.2      
  [7] listenv_0.9.0       fastICA_1.2-3       digest_0.6.31      
 [10] htmltools_0.5.4     fansi_1.0.4         checkmate_2.1.0    
 [13] memoise_2.0.1       cluster_2.1.4       sfsmisc_1.1-14     
 [16] remotes_2.4.2       globals_0.16.2      bdsmatrix_1.3-6    
 [19] prettyunits_1.1.1   jpeg_0.1-10         colorspace_2.1-0   
 [22] skimr_2.1.5         xfun_0.37           callr_3.7.3        
 [25] crayon_1.5.2        jsonlite_1.8.4      survival_3.5-0     
 [28] glue_1.6.2          gtable_0.3.1        sjmisc_2.8.9       
 [31] V8_4.2.2            car_3.1-1           pkgbuild_1.4.0     
 [34] DEoptimR_1.0-11     ggm_2.5             abind_1.4-5        
 [37] scales_1.2.1        rstatix_0.7.1       miniUI_0.1.1.1     
 [40] Rcpp_1.0.10         xtable_1.8-4        htmlTable_2.4.1    
 [43] clue_0.3-64         foreign_0.8-84      Formula_1.2-4      
 [46] stats4_4.2.2        profvis_0.3.7       htmlwidgets_1.6.1  
 [49] RColorBrewer_1.1-3  lavaan_0.6-13       ellipsis_0.3.2     
 [52] farver_2.1.1        urlchecker_1.0.1    pkgconfig_2.0.3    
 [55] nnet_7.3-18         deldir_1.0-6        utf8_1.2.3         
 [58] labeling_0.4.2      tidyselect_1.2.0    rlang_1.0.6        
 [61] reshape2_1.4.4      later_1.3.0         munsell_0.5.0      
 [64] tools_4.2.2         cachem_1.0.6        cli_3.6.0          
 [67] generics_0.1.3      sjlabelled_1.2.0    broom_1.0.3        
 [70] fdrtool_1.2.17      evaluate_0.20       stringr_1.5.0      
 [73] fastmap_1.1.0       yaml_2.3.7          processx_3.8.0     
 [76] knitr_1.42          fs_1.6.0            robustbase_0.95-0  
 [79] glasso_1.11         pbapply_1.7-0       RBGL_1.72.0        
 [82] nlme_3.1-162        mime_0.12           compiler_4.2.2     
 [85] rstudioapi_0.14     curl_5.0.0          png_0.1-8          
 [88] ggsignif_0.6.4      tibble_3.1.8        pbivnorm_0.6.0     
 [91] stringi_1.7.12      ps_1.7.2            lattice_0.20-45    
 [94] Matrix_1.5-3        psych_2.2.9         vctrs_0.5.2        
 [97] pillar_1.8.1        lifecycle_1.0.3     cowplot_1.1.1      
[100] insight_0.19.0      data.table_1.14.6   corpcor_1.6.10     
[103] httpuv_1.6.8        R6_2.5.1            latticeExtra_0.6-30
[106] promises_1.2.0.1    gridExtra_2.3       parallelly_1.34.0  
[109] sessioninfo_1.2.2   codetools_0.2-19    gtools_3.9.4       
[112] pkgload_1.3.2       withr_2.5.0         mnormt_2.1.1       
[115] parallel_4.2.2      rpart_4.1.19        tidyr_1.3.0        
[118] rmarkdown_2.20      carData_3.0-5       shiny_1.7.4        
[121] base64enc_0.1-3     interp_1.1-3